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Abstract. We investigate the equation of state w(z) in a non-parametric form using 
the latest compilations of distance luminosity from SNe la at high z. We combine the 
inverse problem approach with a Monte Carlo to scan the space of priors. On the light 
of the latest high redshift supernova data sets, we reconstruct w(z). A comparison 
between a sample including the latest results at z > 1 and a sample without those 
results show the improvement achieved by observations of very high z supernovae. We 
present the prospects to measure the variation of dark energy density along z by this 
method. 
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1. Introduction 

The acceleration of the rate of expansion of the Universe, first discovered through 
supernovae [TJ [2] is being examined through more extended samples of those distance 
indicators. In the next decade, large samples of supernovae and an increased coverage in 
z of the measurements of baryon acoustic oscillations (BAO) would allow to obtain the 
expansion rate H(z) and the equation of state w(z) [SIS]. Complementary information 
on the matter density and global curvature provided by the CMB measurements from 
Planck [5] and weak lensing surveys [6] would enable to establish a better range of 
allowed regions in w(z). SNe la together with BAO measurements up to very high z 
can be combined to restrict the range of possibilities of the empirical behavior of dark 
energy. By recovering the empirical behavior of w(z) one expects to test whether the 
acceleration of the expansion of the Universe indicates modifications of gravity beyond 
GR, whether it is due to the presence of vacuum energy or it is associated with a light 
scalar field (0[Sl[9l[Tni[ni[l2l[13l[Tl[T5l[l6]). 

The recovery of the function w(z) from a given sample of data has been attempted 
proposing fitting functions or expansion series of w(z) along z in ways to accomodate a 
wide range of dark energy candidates. There has been some debate on the effect that 
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choosing particular models for those functions or truncating the expansion series in z 
might have in deriving possible evolution [T7J [18], [19] . 

Here, we use an approach to obtain w(z) without imposing any constraints on 
the form of the function. This is obtained through a generalized nonlinear inverse 
approach. The inverse approach formulated by Bakus and Gilbert [20] has been widely 
used in geophysics and solar structure physics. In this approach, the mere fact that 
the continuous functional has to be derived from a discrete number of data implies 
the non-uniqueness of the answer. It has also been shown that, even if the data were 
dense and with no uncertainty, there would be more than one solution to many specific 
inverse problems such as the determination of the density structure of the earth from the 
data on the local gravitational field at its surface, and others. This lack of uniqueness 
comes from the way in which the different equations reflect in the observables used. The 
problem of the determination of dark energy faces such degeneracy. In the luminosity 
distance along z from supernovae and other cosmic distance indicators, w(z) enters in an 
integral form, which limits the possibility to access to w(z). In earlier examinations of 
the degeneracy in w(z) obtained through cosmic distance indicators, a range of solutions 
giving the same luminosity distance along z were pointed out [2T] . As more data would 
constrain w(z) at various redshifts the reconstruction should become more successful. 
Here, we examine this using the latest compilations by Wood-Vasey et al (2007) [22] 
and Davis al (2007) [23] which use supernovae gathered by many collaborations and add 
new ones from ESSENCE, SNLS and the Higher-Z Team. 

To compare with a significant body of work which analyses the data using the 
expansion to first order w(z) = w + w a z/(l + z) [2U [25], we will also formulate 
this approach for the case of determination of discrete parameters. We examine from 
current data the possibility of determining at present the values of w(z) and its first 
derivative and compare with previous results. It is known (see for instance [IB] ) that 
this approximation does not allow to recover properly the value of w at high z. And it 
has been argued [19j that this expansion might bias w a towards faster evolution. It is, 
however, a very useful approach to restrict quintessence proposals. Given its widespread 
use, we include here the discrete case with the two parameters w and w a . 

In the following section, we introduce the inverse approach and deduce the equations 
for both the continuous and the discrete case. Then, we show the reconstruction of w(z) 
with current SNela data. Finally, we draw our conclusions. 

2. Inverse problem 

2.1. Non-parametric non-linear inversion 

The inverse problem provides a powerful way to determine the values of functional 
forms from a set of observables. This approach is useful when the information along 
a certain coordinate, in our case information on w(z), emerges in observables coupled 
with information at all other z as it happens with the luminosity distance. Dark energy 
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is here addressed using the non-linear non-parametric inversion. Most frequently, when 
the parameters to be determined are a set of discrete unknowns, the method used is a 
least squares. But the continuous case, where functional forms are to be determined, 
requires a general inverse problem formulation. The inverse method used here is a 
Bayesian approach to this generalization [26J. 

We consider a flat universe with only two dominant constituents (at present): cold 
matter and dark energy. Therefore we characterize the cosmological model by the 
density of matter, Q>m, and by the index w(z) of the dark energy equation of state, 

w{z) = W)- () 

The vector of unknowns M has then a discrete and a continuous component, 

M = fe)- (2) 

On the other hand, the observational data are SNe la magnitudes. We have a finite 
set of N magnitudes, m^, and consider the following theoretical equation, the magnitude- 
redshift relation in a flat universe relating the unknowns to the observational data: 

m th {z, Q M , w(z)) = M + 5 hg[D L (z, Q M , w{z))] - 5 log H + 25, (3) 

where Dl is the Hubble-free luminosity distance 



D L (z,Q M ,w(z)) = c(l + z) - ° (4) 

Jo H{z,{Im,w{z)) 

with 



h(z', n M , w {z)) = h q Jn M (i + z'f + n x ( z >) (5) 



n x (z) = n x ex P ^3 £ ^ / 1 ^f ) ) • 



(6) 



We redefine our data and convert the original SNe magnitudes to dimensionless 
distance coordinates y: 

expio (0* + 51ogifo - 25)/5) [* dz' 

Vi = TT, v = / / ( 7 ) 

+ Jo JQ m (1+z'¥ + Q x (z') 



lnl ° 

*vi = —^—yiVm, (8) 

where \Xi = rrii — M is the distance modulus. With this definition, we deal directly with 
a function w ( z )), the only part which depends on the cosmological model. 

After the corresponding transformations, the observables are now described by a 
vector of iV components, and by a covariance matrix, C y . This method can handle 
correlated measurements, where non-diagonal elements C Viy . are different from zero 
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(observations i and j being correlated). But, at present, those have not been estimated 
for the composite samples of distance indicators. We would then use: 

C v,ij = a iAj ( 9 ) 

Similarly, the unknown vector of parameters is described by its a priori value, Mo, and 
the covariance matrix, C . The function describing w(z) is expected to be smooth, and 
this leads to no null covariance between neighboring points in z for w(z). Thus, the 
covariance matrix C has the form: 

Co=f^ r ° ) (10) 

where a choice is made for the non-null covariance between z and z', C w t z \ w ( z i\. This 
choice is taken to be as general as possible. It would define the smoothness required in 
the solution by setting the correlation length between errors in z and z' (this gives the 
length scale in which the function can fluctuate between redshifts). The amplitude of 
the fluctuation of the function is given by the dispersion a w at z. 
Thus for a Gaussian choice, C w ( z ^ w ( z i) is described as: 

C w{z)Mzl) = al exp I ^^—j > ( n ) 

which means that the variance at z equals and that the correlation length between 
errors is A z . Another possible choice for C w ^ >w ^ z i) is an exponential of the type: 



C w{z)Mz >) = al exp , (12) 




while no difference in the results is found by those different choices. 

This is all the information we have beforehand, and with that we are interested 
in determining the best estimator M for M. The probabilistic approach we will use 
incorporates constraints from priors through the Bayes theorem, i.e, the a posterori 
probability density / poS i(M/D) for the vector M containing the unknown model 
parameters given the observed data D, is linked to the likelihood function L and the 
prior density function for the parameter vector as: 



/ post (M/D) aUP/M) ■ f pri0 r(M) (13) 

The theoretical model described by the operator y th , which connects the model 
parameters M with the predicted data D predicted = y tA (M), is to agree as closely as 
possible with the observed data y. Assuming that both the prior probability and 
the errors in the data are distributed as Gaussian functions, the posterior distribution 
becomes: 
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/ post (M/y) a exp[-- (y - y th (M))* C" 1 (y - y*(M)) 



— (M 

2 v 



M )* Co 1 (M - M ) 



(14) 



where * stands for the adjoint operator. The best estimator for M, M, is the most 
probable value of M, given the set of data y. The condition is reached by minimizing 
the misfit function: 



S = - (y - y th (M)Y C y l (y - y*(M)) + 
~ (M-M )* Co 1 (M-M ), 



(15) 



which is equivalent to maximize the Gaussian density of probability when data and 
parameters are treated in the same way. This Bayesian approach helps to regularize the 
inversion. 

Let us now define the operator G represented by the matrix of partial derivatives 
of the dimensionless distance coordinate, which will simplify subsequent notation. Its 
kernel will be denoted by g as defined in the next equations. 

/ F)„th a n ,th \ 



G 



8Um 


dw(z) 


dyf 


dy\ h 




dw(z) 








dy% 




dw(z) 



(16) 



with 



th 



dVt M 

Qyth 



1 r« il + z'fdz' _ r* , 
= / gn M {z)dz, 
Jo 



HHz>) 



1 f*i 30 x (z')ln(l + z')dz' 

2 Jo 



g w (z')dz'. 



(17) 



dw(z) 2 Jo H 3 {z' 

As shown in Eq.[3] or equivalently Eq.[?l the inverse problem is nonlinear in the 
parameters, thus the solution is reached iteratively in a gradient based search. To 
minimize S in Eq.dH one demands stationarity. For the nonlinear case the solution has 
to be implemented as an iterative procedure where [26] : 



M [fc+1] 

(y ^ 



= M + 

y th (M [k] ) 



Co Gnu (C y + G[fc]Co G 



'[*] 
G [fc] (M 



* \-i 



Mr 



(19) 
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Since we are working in a Hilbert space with vectors containing functional forms, the 
above operator products give rise to scalar products of the functions integrated over 
the domain of those functions. The expressions transform into having the products 
rewritten in terms of the kernels of the operators [27J. 

We will indicate the scalar product by " • " and it is defined as it can be seen from this 
example: 

C w ■ = / C w (z, z')g w {z')dz' (20) 



dw(z) Jo 

The components of the vector of unknowns M, which in our case will be both Qm 
and w(z), are then obtained from: 



N z . 

M [k+1] (z) = M (z) + Y, W m / l C (z,z')g i[k] (z')dz', (21) 

i=i Jo 



where 



TV 

W m = Y.(S[ k \) V m (22) 



V =y+ G (M-Mo)-y^(M) 

V m = yi + r g m {z) (M [k] (z) - M (z)) dz - yl h (z h Sl M , w(z)) (23) 



S = C y + G C G* 

i-Zj r-Zi 

Si,j[k] = {C y )i,j + ^ J o g i[k ](z)C (z,z')g j[k] (z')dzdz' (24) 

In the case of the dark energy equation of state and the matter density the 
expressions reduce to 

2 N 9yt h 

^M[k+l] = Qmo + &n M W i [k] "arT - ( 25 ) 

i=l Oll M [k] 



N 



W[ k+ i](z) = w (z) + 22 w i[k] C w (z,z)g w [ k ](z)dz' (26) 
i=i Jo 

where C w (z,z') = C w ^ >w t z i\(z, z'), Wi[ k ] is given by the product (122]) with: 

f) v th Q th 

V, = Vi + j^(n M - Q Mo ) + -4- • (w - w ) - yf{z h Q M , w(z)) (27) 
o\Im ow(z) 

Qyth Qyth Qyth / Qyth \ 

Sij = SuWi + Q^Cn M ^ + • (C w ■ ^ j (28) 

To test the accuracy of the inversion we use the a posteriori covariance matrix. 
It can be shown (see [28, 27J) that for the linear inverse problem with Gaussian a 



7 



priori probability density function, the a posteriori probability density function is also 
Gaussian with mean Eq.[19]and covariance Eq.[29l Although its value is only exact in 
the linear case it is a good approximation here, since the luminosity distance is quite 
linear on the equation of state w(z) at low redshift. 

Cjyj = ( G* C y 1 G + C 1 ) 1 = Co — Co G* S 1 G Co 

= ( I - C G* S _1 G ) C (29) 

In an explicit form, the standard deviations from this covariance read 



&n u - J Cq m - an 



a 2 



M,)-E^-7^^( 5_1 )M7r-7^- c '- (31) 



\\ w ^> ^ w dw(zY Jl ' J dw(z 



where the symbols with tilde are the a posteriori values, whereas the symbols without 
represent the a priori ones. It must be stressed that the uncertainty in the final w(z) 
does depend on the a priori assumption of the uncertainty In fact, the same w(z) could 
depend on the prior. For avoiding a dependence of the result and its error on the prior 
we use Monte Carlo methods later in the analysis. 

There are other parameters which help to interpret the results. From the form 
of Eq.[29] we see that the operator CoG*S _1 G is related to the obtained resolution. 
This is usually called the resolving kernel K(z,z'). The more this term resembles the 
^-function the smaller the a posteriori covariance function is. In fact, in the linear case, 
the resolving kernel represents how much the results of the inversion differ from the true 
model. It equals to be the filter between the true model and its estimated value [20l [26] . 
In any applied case, it is a low band pass filter which depends on the data available and 
the details requested from the model. In a useful way, it can also be expressed in terms 
of the a priori and the a posteriori covariance matrices: 

K = I - Co 1 . (32) 

This expression will be evaluated numerically to quantify the resolution and information 
generated in the inversion. 

2.2. Discrete parameters 

In the previous section we have obtained the results for a set of a continuous function and 
a discrete parameter, but we can also consider the case of various discrete parameters. 
It was pointed out that a succesful parameterization for modeling a large variety of dark 
energy models is obtained by considering w(z) expanded around the scale factor a. The 
earlier parameterization to first order in z given by w(z) = w + w'z proved unphysical 
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for the CMB data and a poor approach to SN data at z ~ 1. For the case of moderate 
evolution in the equation of state, the most simple (two-parameter) description of w(z) 
so far proposed is [25] El] : 



w(z) = w + w a (l - a) 
where the scale factor a = (1 + z) -1 and w(z) turns out to be: 
w{z) = w + w a - 



l + z 



(33) 



(34) 



We use now this particular form for the function w(z) commonly used to study the 
behaviour of dark energy to solve iteratively for wo and w a : 



w [k+i] =w% + al o W t 



i[k] 



dy 



i=l 

N 



dw, 



w a[k+1] =w° a + a 2 Wa Wi 



dy 



0[k] 



th 



i=l 



where 



dy 



tit 



dw 
dyj h 



dw a [k] 

1 r« 3n x (z')\n(l + z')dz' 

2 Jo mV) ' 

2 Jo HHz') 



The a posteriori variance is for these parameters: 



a 



tt'O 



\ 



> J 



dw " Jl ' J dwo w ° 



a„ 



1 



aw/"- 



The equations for Qm are those of section EH] (Eqs.l25l [17] and [301) . 



Q th Qyth 

"Ml (Q-l\ y l ~2 



(35) 
(36) 

(37) 
(38) 

(39) 
(40) 



3. Determination of w(z) 

In the following, we determine w(z) using the inverse approach described above and 
the SNela by Davis et al. [23j which contains samples of high-z supernovae by various 
collaborations (the Supernova Cosmology Project, the High-Z Team, SNLS, the Higher- 
z team and ESSENCE). 

Explicitly, one obtains the value of Qm and w at a given redshift with the equations 
of section 12.11 (Eqs. [2"5], 1261 [301 and [31]). In all the cases, it is assumed a flat universe, 
where the equations have been deduced, and a good knowledge for the density of matter: 
VL M = 0.27 ± 0.03. 
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Figure 1. (a) Inversion of 1000 data sets generated by bootstrap resampling of the 
latest high redshift compilation in [23j . Initial priors on the equation of state are 
randomly distributed between —3 < w(z) < 1. The solid black line represents the 
mean of the 1000 inversions and the filled regions are the intervals where la and 2a 
of the results lie. (b) The same as (a) for the redshift interval of best reliability. The 
panel shows w(z) (solid line) and the la (dashed shadow) and 2a (pale dashed shadow) 
confidence intervals. 



The a priori model as we will show is arbitrary for w(z), but determines where 
the solution is searched. To avoid the possibility of having converged to a secondary 
minimum just because the prior is too far from the absolute one, we have also carried out 
a Monte Carlo (MC) exploration of the a priori model space. A flat distribution of 1000 
models between —3 < w(z) < 1 has been randomly generated, and then the iterative 
process of inversion has been applied for all of them. We used for this inversion 1000 
data sets obtained by bootstrap resampling of the original one. The MC exploration 
and the bootstrap approach allows to obtain the variance and deviation at every z in 
the non-parametric function w(z). The resulting la is obtained in this way. 

Our results show with a Monte Carlo exploration with priors ranging from w(z) = 
—3 to w(z) = +1 that the obtained solution is stable. 

The solution using the 192 supernovae (60 from ESSENCE, 57 from SNLS and 45 
nearby ones and 30 by Riess et al. 2006) is shown in Figured! The calculation uses 
redshift intervals of Sz = 0.06. Various functional forms for the covariance have been 
tried giving the same results (covariance as in Eq.[12]or the Gaussian Eq.[TT]give similar 
results). This calculation does not restrict fast evolving w(z) as the correlation length 
is kept to allow fast changes of slope (A^ = 0.08). In a situation where the data are 
scarce with very wide priors in the function w(z), the solution can iterate between saddle 
points and local minima as few data do provide a landscape with no strong minima. 
However, in our case, the sample available for w(z) is large enough to allow to find the 
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solution. 

As it has been seen in Figure 1, at low and intermediate redshift the data delimit a 
narrow band of possible solutions. At low and intermediate redshift, where the inversion 
is reliable, there is no sign of any significant evolution. The cosmological constant is 
always within la intervals. 

In Figure 1 the black solid line indicates the evolution of the equation of state, 
whereas the shadowed regions represents the la and 2a intervals. In Figure 2 the 
resolving kernels at z = 0,0.24,0.48,0.84 and 1.20 are shown. It can be seen that 
only the prior is recovered at high redshift, where there is no information coming from 
data and the resolving kernels are almost flat. While Figure 1 show the results for the 
inversion using the [23] sample, very similar results are found if using the sample from 
|22j. The resolving kernels (Figure 2) indicate that the function is generally not well 
resolved at individual redshifts, well beyond the redshift range z ~ 0.5. At high redshift, 
z = 1.2 for example, we observe a very wide and extremely flat K(z, 1.2), meaning that 
this redshift is not resolved at all by the data. The reliability of the inversion peaks in 
the range of z ~ 0.2-0.5 where the information is maximum. This is also found in other 
analyses, though most of the w(z) reconstructions [2U [29], [191 EQ] have been done prior 
to the availability of these combined samples. 

Finally, one can particularize the continuos form of w(z) to a parameterization and 
compare the results with other analyses. Using the equations from Section l2T2l for Qm , w$ 
and w a with priors Om = 0.27±0.03, w = — 1±10 and w a = 0±10 we obtained a positive 
evolution of the equation of state for the full data set: wo = —1.1 ± 0.3, w a = 0.6 ± 1.6 
(when excluding the highest z HST supernovae w = —0.7 ± 0.4, w a = — 2 ± 2). This 
result agrees with the continuous evolutions seen before, as it is expected for a smooth 
and globally monotonous behaviour. For a constant equation of state, both data sets 
favour a cosmological constant within its la intervals with a value of Wq = —1.01 ±0.13 
in the first case and wq = —1.08 ± 0.14 in the second one. 

4. A running cosmological constant as an inverse problem 

We now use the power of the inverse method to detect an evolution of the cosmological 
constant. We interpret the unknown function as a running cosmological constant, but 
it can be interpreted as well as a general function of the dark energy density. 

A running lambda can occur in various different scenarios (see for instance [TTJ [32] and 
references therein). We want to determine a general function AQ A (z) such that 

Q A (z) =n° A + AQ A (z) (41) 

where Q A is the value of Q A (z) for z = 0. We use equations Eq. 0] and [5] to obtain 
the G matrix. We calculate the form of the best tt\(z) by computing dy th /dAQ A (z) in 
iterations until numerical convergence is obtained. 

We tested various degrees of prior knowledge on Qm- We display here the situation 
where we have a very precise knowledge of Qm as it is expected after the Planck mission. 



11 



0.30 



0.25 



0.20 



0.1 5 



0.1 



0.05 



0.00 



1 I.I 1,1 1 1 1 1 1 1 1 1 1 


(a) ; 


: ': h 
'■ : »' \ 
: / 1 \ 

: 1 1 1 ; 
1 1 ] 
1 ' 1 1 : 


- 

■ 


■ 1 1 ' i 1 
: : 1 1 1 j 
:' :i 1 l 1 
■ i I \ ; 
t. 1 i i 
~ "; I: 1 l : 
> : 1 1 1 j 




' : 1 l i : 
'. : 1 i 

SJ : ' \ 1 y v k i 

X ,' V / 1 • • • \ : 








-: \ ' '!*''•■•'' 1 ! \'\ ! 







0.0 



0.5 



1 .0 



1 .5 



0.30 



(b) 



0.25 



0.20 



0.1 5 



0.1 



0.05 



0.00 



0.0 



0.5 



1 .0 



1 .5 



Figure 2. Top panel. Reconstruction of w(z) using the 162 SNe in Ref 22J. The 
different resolving kernels at z = 0,0.24,0.48,0.84,1.20 are shown. The resolving 
kernels at high z show that there is no information to conclude on the evolution of 
the equation of state, there. Bottom panel. The same as before but for the sample 
in Ref. [23] . That allows to see the effect of the HST supernovae on the evolution at 
intermediate redshift. 

If we knew well the matter density, we would see the evolution of Q^(z), as in Figure 3. 

The procedure for the inversion is the same as for the barotropic index of the 
equation of state. The uncertainty on the prior on AQ\(z) is set to a(z) = 0.1. The 
Monte Carlo exploration of the f2A(^)-space is made in the range —0.2 < A^a(^) < 0.2. 
Since we are using the same data sets as in previous sections, we expect the same 
resolution. Therefore we use the same grid resolution. 

For the two samples, VW07 and D07, the result of the inversion is quite similar. 
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Figure 3. Reconstruction of Oa(z) for a flat Universe with Q,m = 0.27. The data 
compilation in VW07 is used in the inversion shown in panel (a). The data compilation 
in D07 is used in the inversion shown in panel (b). The possibilities opened by precise 
measurements of CIm and further data at very high z are illustrated here. 



A constant density cannot be discarded. The reconstruction can not provide further 
information beyond redshift 0.6, but gives the behavior in the range to 0.6. If we relax 
the prior knowledge of Qm to a range of 0.03 we find a similar mean behavior but with 
a wider range of allowed values for Q\. Figure 3 is shown for purposes of indicating the 
enormous value of the complementary information to be provided in the next future. 

These reconstructions can be compared with the behavior of the cosmological 
constant in the running scenarios. Up to now, the obtained smooth evolutions down to 
redshifts lower than 1 are compatible with various scenarios for running cosmological 
constant as no fast changes are expected in the past Gyrs. Scenarios of moderate 
evolution with v ~ 0.1 [32] can be easily discarded with a good knowledge of Qm- 

5. Summary and conclusions 

We introduce here an Inverse Problem approach to determine w(z) as a continuous 
function in a model-independent and non-parametric way. The resulting algorithm 
retrieves w(z) without imposing any constraints on the form of the function. The 
method uses Bayesian information such as the area where this solution is to be found, 
which can be quite unrestrictive. The approach explored here enables to see at which z 
is the maximum information on w(z) or AQ\ in various samples. 

The exploration of dark energy applying this method to the present SNe la sample 
helps to answer the question on whether there is evidence in the evolution of w(z) along 
z. 

With the sets of 162 SNela and 192 SNe la [231 [22] we find w(z) compatible at 
la with a cosmological constant. The larger data set [23] shows improvement in the 
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resolving kernels at high z. Error bars in both cases are big enough to make evident the 
limited knowledge that we have on w(z). However, the bulk of data to come in the next 
decade and the complementary probes would allow to draw the behavior of w(z) and 
AQ\(z) with much higher precision. This method can be applied to the determination 
of w(z) with both SNe la and BAO measurements. 
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